function [x,normk,it,nrm] = solverNT_mc(in,xo)
	% this function solves the coarse mesh problem via what is essentially
	% a classical newton method.  Here, we can compute nearly completely the
	% Jacobian given the definition of the nonlinear function.  Only values
	% which are explict functions of keff need to be approximated via finite
	% differences.  For reference, the nonlinear function is:
	%  FF = [ R*M-x(end)*eye(length(R(:,1)))  zeros(length(R(:,1)),1) 
	%         F*M                            -sum(A*M*x(1:end-2))-L          ];  
	% z1 = FF*(x(1:end-1)); 
	% z2 = 0.5 - 0.5*sum(x(1:end-2));

	disp('-------------------------------------------')
	disp('--      beginning newton iterations      --')
	disp('-------------------------------------------')

	normk = 0;
	x           = xo;%;[xo' 1.0]'; % add in current eigenvalue
	x(end)=1;
	%x(end-1)=0.98988;
	xx(:,1)     = x;
	[z,R,F,A,L] = resid2D(x,in); % initial residual
	itmx        = 20; 
	it          = 1;
	nrm(1)      = norm(z);

	while ( nrm(it) > in.tol(1)*nrm(1)+in.tol(2) && it < itmx )
	    disp(['  iteration number: ', num2str(it), ' norm(x) = ', num2str(nrm(it))])
	   % fp          = jacob2D(x,in,z); % FD appx (for verification)
	    fp          = jacob2Da(x,in,z,R,F,A,L); % hand coded
	    s           = -fp\z;
	    x           = x + s;
	    it          = it+1; 
	    xx(:,it)    = x; 
	    [z,R,F,A,L] = resid2D(x,in);
	    nrm(it)     = norm(z);
	end

end % solverNT

function [z,R,F,A,L] = resid2D(x,in)
	%% function z = resid2D(x,input)
	% This function computes the nonlinear residual
	% Here, x(1:end-2) corresponds to the legendre coefficients
	%       x(end-1) is keff
	%       x(end)   is the current eigenvalue (lambda)
	% NOTE: I return the most recent response functions
	M     = in.M;
	order = in.order;
	numel = in.numel;
	numg  = in.numg;

	[R,F,A,L] = combmesh_center_2d_resp(in,x);
        %disp('-------')
        %size(L)
        %size(A)
        %size(L*M*x(1:end-2))
        %size(x(1:end-2))
        %disp('-------')
	loss = (A+L)*M*x(1:end-2);

	%compute the nonlinear function, FF(x)*x = 0
	FF = [ R*M-x(end)*eye(length(R(:,1)))  zeros(length(R(:,1)),1)
	       F*M                            -loss  ];
	z1 = FF*(x(1:end-1));
	z2 = 0.5 - 0.5*x(1:end-2)'*x(1:end-2);
	z  = [z1' z2']';
end % resid2D

function fp = jacob2Da(x,in,z,R,F,A,L)
	%  x  = current unknowns vector
	%  in = problem data
	%  z  = last computed residual
	M = in.M;
	n   = length(x);
	eps = 1.e-8;
	fp  = zeros(n);

	%[R,F,A,L] = combmesh_center_2d_resp(in,x);

	fp(1:end-2,1:end-2) = R*M - x(end)*eye(n-2);
	fp(n-1,1:n-2)       = (F-x(end-1)*(A+L))*M;
	fp(n,1:n-2)         = -x(1:end-2)';
	fp(1:n-2,n)         = -x(1:end-2);

	% compute the perturbed response functions
	x(end-1) = x(end-1)+eps;

	[R2,F2,A2,L2] = combmesh_center_2d_resp(in,x);

	Rk = (R2-R)/eps;
	Fk = (F2-F)/eps;
	Ak = (A2-A)/eps;
	Lk = (L2-L)/eps;
	fp(n-1,n-1)     = (Fk-x(end-1)*(Ak+Lk)-(A+L))*M*x(1:end-2);
	fp(1:n-2,n-1)   = Rk*M*x(1:end-2);
end % jacob2D

function fp = jacob2D(x,in,f0)
	% fp = jacob(f,x,fo)
	%    this function produces a crude forward finite difference
	%    approximation for the function f(x) about x; fo is the initial
	%    evaluation of f
	n   = length(x);
	%f0  = feval(f,x);
	eps = 1e-8;
	fp  = zeros(n);
	for j=1:n
	    v           = zeros(n,1); 
	    v(j,1)      = 1.0;
	    xx          = x+v*eps;
	    f1          = resid2D(xx,in);
	    fp(:,j)     = (f1 - f0)/eps;
	end
end

function [R,F,A,L] = combmesh_center_2d_resp(in,x)
M = in.M;
order = in.order;
numel = in.numel;
numg = in.numg;
% initialize response matrices
R = zeros( numel*numg*(order+1)*4 );
F = zeros( 1, numel*numg*(order+1)*4 );
A = zeros( 1, numel*numg*(order+1)*4 );
L = zeros( 4, numel*numg*(order+1)*4 );
lenR = numg*(order+1)*4;
lenM = length(M);
% build mm
%mm = zeros(length(M),1);
for i = 1:length(M(:,1))
    if ( sum(M(:,i)) + sum(M(i,:)) )==0
        mm(i)=1;
    else
        mm(i)=0;
    end
end
mm = mm(1:numg*(order+1):end);
mm = spdiags( mm', 0, length(mm), length(mm));

% -------------------------------------------------------------------------
% Generate Response Matrices

[Rb Fb Ab Lb]  = mesh_center_2d_resp(x(end-1),in);

% this gets the info for each element type; now, need to place
I = length(in.elements(:,1));
J = length(in.elements(1,:));
k = 0;

Rlist = cell(1,in.numel);
Llist = cell(1,in.numel);
for j = 1:J      % rows
    for i = 1:I  % cols
        
        if in.elements(i,j) > 0
            % I am the kth element
            k = k + 1;
            type = in.elements(i,j);
        else
            % else I'm just void
            continue
        end
        tmp(1:lenR,1:lenR)=Rb(type,:,:);
        Rlist{k} = sparse(tmp);
        %R( ((k-1)*lenR+1):(k*lenR),((k-1)*lenR+1):(k*lenR) ) = Rb(type,:,:);
        F(1, ((k-1)*lenR+1):(k*lenR) ) = Fb(type,:);
        A(1, ((k-1)*lenR+1):(k*lenR) ) = Ab(type,:);

        myL(:,:) = Lb(type,:,:);
        my1 = myL(1,:);
        my2 = myL(2,:);
        my3 = myL(3,:);
        my4 = myL(4,:);
        
        % for an element, we have
        Llist{k}  = sparse([my1 my2 my3 my4
                            my2 my1 my4 my3
                            my4 my3 my1 my2
                            my3 my4 my2 my1]);
        % build the L matrix
       % LK( (k-1)*4+1:(k-1)*4+4, (k-1)*4*tmp+1:k*4*tmp ) = Ltmp;
      %  LKlist{k} = Ltmp;
    end
end

% Build R
R   = blkdiag(Rlist{:});
LK  = blkdiag(Llist{:});
tmp = numg*(order+1);


% NOTE: 3.16.2010 -- I added the M because the input is actually J+
%  and LK is meant to couple with J-
L = mm*LK;
%size(L)
%figure(2)
%spy(L)
%print -djpeg90
%disp('old leakage')
%sum( mm*( (M*x(1:end-2))'*LK')' )
%disp('new leakage')
%sum(L)*M*x(1:end-2)
L=sum(L);
end % combmesh_center_2d_resp





